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Abstract 


This  project  analyzes  the  effect  of  transverse  shear  deformation  upon  the  aeroelastic  response 
of  composite  wings  in  high  speed  flow  regimes.  Previously,  models  have  been  developed  to 
predict  the  aeroelastic  characteristics  of  classical  materials  in  high  speed  flow.  However, 
these  studies  ignored  transverse  shear  by  assuming  an  infinite  modulus  of  rigidity.  This 
assumption  underestimates  transverse  flexibility  by  ignoring  the  transfer  of  loads  through 
the  wing  thickness.  By  assuming  a  finite  modulus  of  rigidity  and  redeveloping  the  governing 
equations,  this  model  would  more  accurately  predict  the  aeroelastic  response  of  composite 
wings.  This  present  analysis  concerns  mainly  the  determination  of  aeroelastic  trends  vice 
more  detailed  solutions.  Thus,  linearized  flow  theory  is  used.  Upon  conclusion,  this  study 
gives  results  for  divergence  speed  and  flutter  speeds,  as  well  as  their  mode  shapes. 
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Chapter  1 
Background 


Since  the  days  of  the  Wright  Flyer,  the  aviation  industry  has  constantly  searched  for  lighter, 
stronger,  and  more  durable  materials.  In  the  past  century  this  pursuit  has  taken  aeronautical 
engineers  from  cloth  and  wood  to  paper  thin  steel  to  honeycombed  aluminum  structures. 
The  next  logical  step  in  this  evolution  is  composite  materials.  However,  due  to  the  differences 
in  the  properties  of  composites  and  the  so-called  classical  materials,  designers  need  guidance 
to  predict  the  characteristics  of  composite  aeroelastic  structures. 

The  purpose  of  this  project  is  to  develop  a  model,  which  will  predict  the  impact  of  using 
composite  materials  in  aircraft  wings.  In  particular,  it  will  look  at  the  aeroelastic  effects 
caused  by  transverse  shear  deformation. 

1.1  A  Brief  Overview  of  Aeroelasticity 

Often,  the  complexities  of  design  force  engineers  to  break  down  their  work  into  simpler 
components.  In  aeronautical  engineering,  the  two  most  vital  of  these  components  are  the 
aerodynamic  and  structural  characteristics  of  the  aircraft.  These  basic  pieces  of  the  puzzle 
answer  the  questions:  a)  will  the  aircraft  fly?  and  b)  will  it  stay  in  one  piece?  In  a  perfect 
undergraduate  level  world,  answering  these  two  questions  would  be  enough.  However,  in 
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the  real  world  things  are  not  as  simple.  In  truth  the  aerodynamic  and  structural  forces 
on  an  aircraft  are  dependent  upon  each  other.  This  being  said,  the  next  logical  question 
becomes,  “How  do  the  two  affect  each  other?”  In  response  to  this,  engineers  began  the  study 
of  aeroelasticity. 

Aeroelasticity  can  be  dehned  as  the  study  of  the  interaction  between  the  aerodynamic, 
inertial,  and  structural  forces  acting  on  an  object. Taking  principles  from  the  helds  of 
aerospace  and  mechanical  engineering,  the  aeroelastician  is  primarily  concerned  with  de¬ 
termining  the  effects  of  placing  an  aerodynamic  load  on  a  structure.  Although  usually 
associated  with  aircraft  design,  aeroelasticians  are  used  in  a  variety  of  fields.  For  instance, 
the  most  famous  of  all  aeroelastic  failures  occurred  on  the  Tacoma  Narrows  Bridge,  which 
can  be  seen  in  Fig.  1.1.  In  1940,  the  aerodynamic  load  caused  by  wind  blowing  through  the 


Figure  1.1:  Torsional  oscillation  of  Tacoma  Narrows  Bridge.^ 


CHAPTER  1.  BACKGROUND 


11 


valley  caused  this  suspension  bridge,  constructed  of  concrete  and  steel,  to  twist,  bend,  and 
eventually  collapse,  as  if  it  were  made  of  rubber.  The  failure  of  the  Tacoma  Narrows  Bridge 
can  be  attributed  to  what  is  known  as  flutter  instability.  Although  not  all  aeroelastic 
events  are  as  visually  dramatic,  they  do  occur,  and  ignoring  them  can  potentially  lead  to 
catastrophic  results. 

The  study  of  aeroelasticity  can  be  broken  down  into  two  major  branches:  static  and 
dynamic.  The  case  of  static  aeroelasticity  is  concerned  with  systems  in  equilibrium.  Once 
an  aerodynamic  load  is  placed  on  an  aircraft,  its  structures  will  deform,  redistributing  the 
load.  This  can  lead  to  one  of  two  possibilities.  The  first  is  simply  a  new  state  of  equilibrium 
in  which  the  aerodynamic  and  structural  characteristics  of  the  wing  are  slightly  changed  for 
better  or  worse.  The  second,  and  less  appealing,  case  is  that  the  redistributed  loads  will 
escalate  until  the  wing  fails.  This  is  known  as  static  divergence. 

Dynamic  aeroelasticity  is  concerned  with  time  dependent  instabilities.  These  can  be 
either  transient  or  oscillatory  depending  upon  the  nature  of  the  response.  The  most  promi¬ 
nent  of  these  cases  is  flutter,  which  simply  represents  a  harmonically  oscillating  wing  and 
constitutes  the  stability  boundary  between  damped  and  undamped  oscillations.  Aeroelastic 
methods  can  be  used  to  predict  whether  the  system  will  eventually  stabilize  or  diverge,  as 
in  the  case  of  the  Tacoma  Narrows  Bridge.^ 

1.2  Composites  vs.  Classical  Materials 

In  recent  years,  a  continual  improvement  in  composite  materials  has  given  engineers  a  new 
level  of  freedom  in  design.  These  new  materials  allow  for  lightweight,  high  performance 
structures,  which  could  not  have  been  constructed  from  metals.  Metals  are  known  as  the 
“classical”  materials.  However,  with  this  change  from  metallic  to  composite  structures,  the 
classic  structural  models  must  be  reexamined  to  determine  if  they  can  be  accurately  applied 
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to  these  new-age  materials. 

The  ability  to  tailor  materials  to  increase  their  functionality  truly  lends  itself  to  the 
held  of  aircraft  and  spacecraft  design.  In  the  quest  to  go  faster  and  hy  higher,  weight 
and  structural  stability  are  vital  components.  Over  the  past  two  decades,  the  effects  of 
these  materials  on  the  aeroelastic  behavior  of  aircraft  wings  have  been  examined.  Although 
composite  technology  is  allowing  for  the  creation  of  wing  structures  of  enhanced  efficiency, 
the  incorporation  of  the  new  technology  is  forcing  aeroelasticians  to  look  back  upon  the  old 
models. 

In  many  cases,  the  classical  models  developed  over  the  past  sixty  years  need  only  to  be 
extended,  as  the  composite  materials’  behavior  closely  resembles  that  of  the  metals.  However, 
composites  tend  to  differ  greatly  from  metals  in  the  case  of  transverse  shear  flexibility. 
Shear  effects  are  created  from  the  transfer  of  loads  through  the  thickness  of  a  structure. 
The  greater  the  shear  rigidity,  the  less  a  material  will  transversely  deform  under  a  given 
load.  Early  on  in  the  derivation  of  the  classical  model,  the  assumption  of  inhnite  rigidity  in 
transverse  shear  was  made.  This  assumption,  referred  to  as  Kirchhoff’s  hypothesis,  can  be 
made  and  justihed  in  metallic  structures.  Conversely,  composite  materials  have  been  shown 
to  have  a  much  lower  modulus  of  transverse  shear  rigidity. *^4, 10)  In  some  cases  the  effects 
of  transverse  shear  deformation  have  been  shown  to  affect  the  static  aeroelastic  response  of 
a  composite  wing  by  almost  hfty  percent. This  is  obviously  not  negligible.  Due  to  the 
large  impact  of  transverse  shear  upon  composite  materials,  the  classical  model  cannot  be  used 
or  even  extended  to  accurately  predict  the  aeroelastic  characteristics  of  composite  wings,  but 
rather  a  new  model  must  be  created  taking  into  account  a  hnite  rigidity  in  transverse  shear 
deformation. 


Chapter  2 
Procedure 


In  order  to  accurately  develop  a  mathematical  model  representative  of  an  aircraft  wing  in 
high-speed  flow,  a  four-step  process  is  required.  First,  the  equations  of  motion  must  be 
developed  for  the  system.  Second,  the  structural  mechanics  of  the  wing  must  be  inserted 
into  these  equations.  Next,  the  aerodynamic  loads  must  be  included  to  create  a  system  of 
governing  equations.  Finally,  this  system  must  be  solved  to  determine  the  critical  aeroelastic 
eigenvalues  and  mode  shapes  of  the  wing.  In  order  to  accomplish  this  last  step,  the  differential 
equation  solver  in  MATHEMATIC  A  was  used. 

2.1  Equations  of  Motion 

The  first  step  in  analyzing  any  physical  system  is  to  choose  a  set  of  axes,  which  can  accurately 
represent  that  system.  Fig.  2.1  illustrates  the  geometric  model  of  a  generic  swept  wing.  The 
wing  geometry  is  described  by  a  Cartesian  coordinate  system  with  xi  set  as  the  effective 
wing  root,  X2  set  normal  to  Xi  as  the  reference  axis  along  the  wingspan,  and  x^  describing 
the  normal  direction  to  the  wing  surface.  This  three  dimensional  system  will  be  used  to 
describe  the  wing  displacement  under  static  and  dynamic  loads  resulting  from  airflow  over 
the  wing.  A  is  used  to  describe  the  angle  of  sweep  with  positive  angles  associated  with  swept 
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Figure  2.1:  Geometry  of  a  generic 

back  wings  and  negative  angles  describing  a  forward  sweep.  This  becomes  significant  as  Ref. 
[5]  shows  that  even  in  low  speed  flow,  the  angle  of  sweep  has  a  dramatic  effect  on  transverse 
shear.  Due  to  the  fact  that  span  is  mnch  greater  than  chord  and  thickness  in  most  common 
wings,  the  three  dimensional  wing  will  be  rednced  to  a  one-dimensional  system  with  only  a 
spanwise  variation  of  properties. 

Having  dehned  the  coordinate  system  to  be  used,  the  equations  of  motion  can  now  be 
developed.  Assnming  hrst-order  transverse  shear  deformation,  that  is  the  displacement  held 
varying  linearly  throngh  the  thickness,  the  three-dimensional  time  dependent  displacement 
eqnations  are  given  by: 
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Ui(xi,  X2,  x^;  t)  =  Ui{xi,  X2]  t)  +  x^.'lplixi,  X2]  t)  (2.1) 

U2{xi,  X2,  x^;  t)  =  U2{xi,  X2]  t)  +  X2,'4)2{Xi,  X2,  t)  (2.2) 

t/3(a;i,  X2,  X3]  t)  =  u:i{xi,X2]  t)  (2.3) 


In  Equations  (2.1-2),  the  first  terms  on  the  right-hand  side  represent  the  displacement 
components  in  the  reference  plane  (0:3  =  0),  while  the  second  terms  represent  the  displace¬ 
ment  off  the  reference  plane.  '^1  is  twist  about  the  X2  axis,  and  1^2  is  twist  about  the  xi 
axis.  Equation  (2.3)  shows  that  the  normal  displacement  of  any  point  in  the  wing  structure 
is  assumed  to  be  the  same  as  a  point  in  the  reference  plane.  Thus,  there  is  a  constant  wing 
thickness  during  deformation. 

To  further  simplify  the  system,  it  is  assumed  that  chordwise  rigidity  exists.  Thus,  ui  0, 
U2  — ^  U2{x2]  t),  and  '^1  — >  'ipi{x2',  t),  which  is  dehned  as  9{x2]  t),  twist  about  the  pitching  axis. 

The  next  step  is  to  rewrite  this  three-dimensional  displacement  field  as  a  one-dimensional 
system.  This  is  a  valid  assumption  as  long  as  the  wing’s  span  is  much  larger  than  its  chord, 
as  it  resembles  a  beamlike  structure;  thus,  the  larger  the  aspect  ratio,  the  more  accurate  the 
model.  The  rotation  about  the  Xi  axis,  'ip2{xi,  X2;t),  is  modeled  as  a  linear  function  in  Xi.  'ip2 
is  the  rotation  of  the  reference  axis  about  Xi,  while  1^2  is  the  rotation  of  the  wing  elements 
off  the  reference  axis.  This  can  be  written  as 

^|J2{Xl,  X2;  t)  =  ^2(2^2;  t)  +  XI'^2{X2-,  t)  (2.4) 

The  normal  displacement  of  the  reference  plane  can  be  written  as 

U‘i{xi,X2]  t)  =  h{x2]  t)  -  {xi  -  xo)9{x2;  t)  (2.5) 


where  h{x2',t)  is  the  vertical  displacement  due  to  plunging  as  measured  at  the  elastic  axis. 
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Inserting  Equations  (2.4-5)  into  Equations  (2.1-3)  results  in  the  one-dimensional  dis¬ 
placement  components.  These  equations  are  now  only  a  function  of  time  and  the  spanwise 
coordinate,  and  can  be  written  as 


Ui  =  Xs9 

(2,6) 

U2+X3{p2+Xi^2) 

(2.7) 

=  h  —  {xi  —  xo)6 

(2.8) 

2.2  Hamilton’s  Principle 

Although  seemingly  complex,  the  entire  method  used  to  hnd  the  wing  modes  is  based  on  the 
principle  of  conservation  of  energy.  Basically,  the  kinetic  energy  generated  by  the  unsteady 
aerodynamic  loads  is  transfered  to  the  wing  creating  potential  energy  in  the  form  of  a  strain. 
The  structural  analysis  then  determines  how  the  wing  will  react  under  these  loads,  either 
deforming  into  a  new  state  of  equilibrium  or  failing.  The  equations  of  motion  result  from 
the  application  of  Hamilton’s  variational  principle.  Minimizing  the  function 

dJ  =  0 

dJ=  /*'{-  [  aij5Uijdip+  f  p{Hi  -  Ui)mid^  +  f  aiSUidn}  dt  (2.9) 

\J  J  ip  \J  ip  \J 

The  hrst  term  on  the  right-hand  side  represents  the  strain  energy  present  inside  the 
structure  of  the  wing,  aij  corresponds  to  the  internal  stresses,  and  (p  represents  the  volume 
of  the  wing.  The  second  term  represents  the  body  forces  and  kinetic  energy  of  the  system, 
with  p  being  the  density  of  the  structure.  The  final  term  denotes  the  energy  resulting  from 
the  surface  stresses,  ps-  is  the  wing’s  surface  area. 
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Hamilton’s  Principle  is  an  application  of  variational  calculus  and  states  that  a  physical 
system  will  go  from  one  state  to  the  next  through  the  lowest  possible  change  in  energy.  By 
integrating  over  time  and  setting  the  change  in  energy  of  the  system  equal  to  zero,  the  state 
of  deformation  at  ti  can  be  determined.  Because  the  variational,  6U  does  not  necessarily 
have  to  equal  zero,  its  coefficients  do.  Thus,  from  inserting  the  structural  properties  and 
aerodynamic  forces  into  Equation  (2.9),  the  displacements  at  a  given  state  can  be  found. 

2.3  Structural  Analysis 

It  is  now  necessary  to  insert  the  structural  properties  of  the  wing.  However,  due  to  the  fact 
that  the  displacements  are  still  general  components  at  this  point.  Equation  (2.9)  must  be 
rewritten  in  terms  of  strain  in  order  to  simplify  the  system.  The  following  relations  are  used 
to  accomplish  this.*^’^’^"^’^^ 


eii  =  t/i,i  =  0  (2.10) 

€22  =  ^2,2  =  U2  +  X3(i^2  +  (2-11) 

€33  =  ^3,3  =  0  (2.12) 
7l2  =  f^l,2  +  f^2,l  =  X3(0'  +  1P2)  (2.13) 

7i3  =  Ui^3  +  =  0  (2-14) 

723  =  h^2,3  +  f^3,2  =  '*(’2  +  Xitl)2  +  h'  +  [XqO)'  —  Xi9'  (2-15) 


Here  the  subscripts  {i,j)  indicate  differentiation  of  the  component  with  respect  to  the 
jth  variable.  Of  particular  interest  to  this  model  is  Equation  (2.15).  This  is  where  transverse 
shear  deformation  can  be  accounted  for.  Had  723  =  0,  Kirchoff’s  Hypothesis  would  be 
satished  as  713  also  equals  zero.  Recall  that  Kirchoff’s  hypothesis  is  the  assumption  made 
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for  classical  materials  in  which  transverse  shear  rigidity  is  inhnite,  causing  transverse  shear 
strains  to  vanish. 

After  the  strain  relations  are  inserted,  Equation  (2.9)  becomes  seemingly  more  complex. 
In  order  to  more  easily  deal  with  it,  the  following  relations  are  dehned 


J  aijX^^x'^  dA 

(2.16) 

[  pHiX^x^  dA 

A 

(2.17) 

f  px^x^  dA 

J  A 

(2.18) 

Equations  (2.16-8)  represent  the  generalized  stress  couples,  body  forces,  and  mass  respec¬ 
tively.  From  these,  the  most  general  equations  of  motion  can  be  derived.  These  are® 


Su2  :  =  0  (2.19) 

6^2  ■  =  0  (2.20) 
5^2  :  +  ^23’°^  -  =  0  (2.21) 

6h  :  -L-  =  0  (2.22) 

69  :  (J(°’2)  +  /(2,o)  _ 

-  xoTS’^^'  -T-  rf =  0  (2.23) 


The  natural  boundary  conditions  also  result  from  Hamilton’s  variational  principle,  and 
are  given  by 

Root  conditions  {x2  =  0): 


U2  =  il)2  =  '4^2  =  h  =  6  =  ^ 


(2.24) 
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Tip  conditions  {x2  =  /): 


.rTi(0,0)  .rTi(0,0) 

-^22  —  -^22 

(2.25) 

t“'(07)  _ ^(07) 

-‘-22  —  -‘22 

(2.26) 

7  22  —  -‘22 

(2.27) 

r^(0,0)  7^(050) 

-^23  ~  23 

(2.28) 

_  ^(07)  ^(1:0) 

7  23  —  742  7  23 

(2.29) 

The  system,  Equations  (2.19-29),  represents  the  most  general  one-dimensional  system. 
The  only  assumptions  made  up  to  this  point  have  been  chordwise  rigidity  and  constant 
thickness.  Now,  the  assumption  that  the  wing  is  composed  of  a  single  layer  of  composite  is 
made.  This  assumption  eases  the  process  of  determining  the  stress  resultants,  and  if  deemed 
necessary,  future  study  can  be  undertaken  on  wings  with  more  than  one  layer. 

The  next  step  is  to  begin  to  dehne  the  material  properties  of  the  system  in  order  to 
analyze  a  more  specihc  wing.  By  now  inserting  the  constitutive  equations,  the  generalized 
stress  couples  can  be  put  in  terms  of  the  displacement  components.  The  three-dimensional 
constitutive  equations  are: 


/(^ll  \  /  Qll  Qi2  Qi6  \  /  \ 

I  <722  1  =  I  Q21  Q22  Q26  1  I  1  (2.30) 

\<7i2/  VQei  Q62  Qm  J  \712/ 


/  <7i3  A  _  /  C45  C45  A  f  713  a 
\0'23  J  V  C54  C55  J  V  723  / 


(2.31) 


Introduction  of  the  three-dimensional  constitutive  equations  into  the  one-dimensional 


stress  tensors  can  be  done  through  the  following  relations: 


[-^ij  (7:2) )  (^ij  (7:2) )  Ojj 


Aij[l,xi,xl]  dxi 


(2.32) 
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[Bij(x2),bij(x2),bij]  = 

J  Bij[l,xi,xj]  dxi 

(2.33) 

\Bij  (3^2 ) )  dij  (^X2 ) ,  dij  ] 

J  Dij[l,xi,xj]  dxi 

(2.34) 

where 


[Aij{xi,X2),Bij{xi,X2),Cij{xi,X2)]  =  f  Qij[l,X3,xl]  dxs  (2.35) 

J  0 

Equations  (2.30-1)  represent  any  material  containing  monoclinic  symmetry,  i.e.  sym¬ 
metry  with  respect  to  the  vertical  coordinate.  From  these  principles  comes  higher  level 
aeroelastic  theory,  such  as  structural  tailoring.  However,  further  rehnement  of  the  structural 
model  can  be  used  as  the  current  emphasis  is  being  placed  on  general  trends,  not  exact  solu¬ 
tions.  With  this  in  mind,  further  assumptions  can  be  made.  In  order  to  reduce  the  equations 
to  a  more  manageable  size,  it  is  assumed  that  the  structural  properties  are  constant  in  the 
spanwise  direction  and  also  that  the  reference  axis  is  along  the  wing  mid-chord.  Now  the 
generalized  stress  resultants  can  be  written  as: 


=  A22u'^ 

(2.36) 

Tif)  =  ^22^2  +  ^26^2  +  ^26^2 

(2.37) 

^23’°^  =  ^55'*/’2  +  A^bh'2  +  ^55(3^0^)2 

(2.38) 

=  d22^'^ 

(2.39) 

^12’  ^  =  B)q2'^2  +  DqqO'  +  DQQlp2 

(2.40) 

^2,3’*^^  =  h55'?/’2  ~  d^^d' 

(2.41) 

The  hnal  step  in  developing  the  wing’s  structural  model  is  to  replace  the  one-dimensional 
stiffness  quantities  with  the  material  properties.  This  can  be  done  in  terms  of  three  proper¬ 
ties:  Young’s  modulus,  the  modulus  of  rigidity,  and  Poisson’s  ratio.  Young’s  modulus,  E,  is 
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a  measure  of  the  ratio  of  stress-to-strain  in  a  material  and  gives  a  comparison  between  the 
stiffnesses  of  various  materials.  The  modulus  of  rigidity,  G  relates  shearing  stress-to-shearing 
strain  in  much  the  same  manner.  Poisson’s  ratio,  z/,  relates  the  axial  strain  to  the  lateral 
strain.  As  can  be  guessed,  all  three  properties  are  related.  However,  since  each  measures  a 
slightly  different  property,  it  is  necessary  to  include  all  three  in  the  material  analysis.  This 
being  said,  the  one-dimensional  stiffness  quantities  can  now  be  written  as^^’^^^ 


A55  —  tcG 


(2.42) 


055  =  t 


doo  — 


144(1  -  z/2) 
D2(\  =  0 


^  24(1  + o'* 


(2.43) 

(2.44) 

(2.45) 

(2.46) 

(2.47) 


where  t  is  the  wing  thickness,  and  c  is  the  chord  length. 

The  last  set  of  assumptions  to  be  made  concern  the  body  forces  and  inertial  terms. 
Because  the  wing  deformation  is  measured  due  to  the  aerodynamic  loads  and  not  the  weight 
of  the  wing,  the  body  forces  will  be  ignored.  Doing  this  signihcantly  reduces  the  complexity 
of  the  governing  equations  and  boundary  conditions.  Also,  at  this  point  the  system  can  be 
reduced  to  four  governing  equations  and  eight  boundary  conditions  by  substitution. 


With  all  the  structural  pieces  in  place,  the  governing  equations  can  now  be  written  in 


terms  of  the  displacement  components,  the  structural  properties,  and  the  aerodynamic  forces. 
The  governing  system  of  aeroelastic  equations  is  given 
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d'i/j2  ■ 


E 


5^2  : 


E 


12(1  -  z/2) 

E 


ct^'ip2  +  tcG'ip2  +  tcGh'  +  XQtcGO'  =  0 


144(1  -  z/2) 


24(1  +  u) 


3  G  3 
12 


t)ff+[ 


E 


24(1  +  u) 


3  G  3 

cep-cH 


5h  :  tcG'ip2  +  tcGh'  +  xotcGd"  +  L  =  0 

+  xltcG)e”  +  -  PcH 

+XotcG'4!2  +  XotcGh"  +  T  =  0 


with  the  boundary  conditions: 
X2  =  0: 

X2  =  /: 


E 


'ilj2  =  'ilJ2  =  h  =  9  =  0 


=  ^2  =  0 

ij2  +  h'  +  xo9'  =  0 
E 


^24(1 +  z/) 


ct^  -  ^Gt]'ilj2  + 


G 


^24(l  +  .f''  +  l2"'')^'  =  ° 


(2.48) 

0  (2.49) 

(2.50) 

(2.51) 


(2.52) 

(2.53) 

(2.54) 

(2.55) 


After  a  good  deal  of  math  and  material  science,  the  structural  model  for  the  wing  is  now 
complete,  but  this  is  only  the  hrst  half  of  the  development.  The  next  step  in  the  process 
is  to  blend  the  aerodynamics  into  this  model.  This  blend  of  disciplines  is  what  makes  the 
study  of  aeroelasticity  difficult,  and  at  the  same  time,  interesting. 
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2.4  Static  Aerodynamic  Forces 

In  order  to  complete  the  governing  eqnations,  the  aerodynamic  load  and  torque  per  unit  span 
must  be  inserted  into  the  governing  equations.  Due  to  the  fact  that  structural  instability 
is  the  amplihcation  of  small  deformations  in  the  wing,  the  aeroelastic  response  to  large 
deformations  represents  a  post-stability  analysis,  which  is  unnecessary.  This  allows  for  the 
use  of  linear  aerodynamic  theory.  Although  not  exact,  the  linear  equations  for  load  and 
torque  provide  enough  accuracy  for  this  study. 

In  aerodynamic  theory  there  are  two  very  different  types  of  flow:  steady  and  unsteady. 
Steady  flow  can  be  assumed  along  constant  flight  paths  through  laminar  fluids,  and  its 
properties  are  independent  of  time.®  Obviously  from  the  mathematical  standpoint,  this  is 
the  more  desirable  type  of  flow.  From  these  equations,  the  static  aeroelastic  equations  can 
be  determined.  Solving  this  system  allows  for  calculation  of  the  flow  velocity  at  which  the 
static  instability,  known  as  wing  divergence,  will  occur,  as  well  as  the  wing  mode  shapes 
obtained  for  sub-critical  flow. 

In  the  case  of  steady  flow,  lift  per  unit  span  can  be  obtained  by 

L{x2)  =  qnac{ao  +  9  —  h'tanA)  (2.56) 

Equation  (2.56)  shows  that  much  of  the  lift  is  determined  by  known  properties  of  the 
system,  such  as  dynamic  pressure,  lift  curve  slope,  chord  length,  and  angle  of  attack.  How¬ 
ever,  notice  that  9  and  h'  come  into  this  equation.  These  terms  represent  the  aeroelastic 
interaction  between  the  aerodynamics  and  the  structural  properties.  The  terms  in  parenthe¬ 
sis  together  are  known  as  the  effective  angle  of  attack.  This  term  accounts  for  the  angle  of 
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attack  at  the  wing  root  as  well  as  any  bending  and  twisting,  which  may  have  occured  along 
the  span. 

For  incompressible  flow,  the  lift  curve  slope  is  only  a  function  of  the  geometric  properties 
of  the  wing.  In  compressible  flow,  a  also  becomes  a  function  of  Mach  number.  By  varying  the 
Mach  number,  and  thus  the  lift  curve  slope,  the  effects  of  high  speed  flow  can  be  evaluated. 
Likewise,  the  torque  per  unit  span  can  be  found  from: 

T(x2)  =  qnace(ao  +  9  -  h'tank)  +  qnC^Cmac  (2-57) 

Definitions  of  the  aerodynamic  properties  can  be  found  in  A.l. 

Replacing  L  and  T  in  Equations  (2.48-51)  with  Equations  (2.56-7)  gives  the  complete 
static  aeroelastic  system  of  equations  accounting  for  both  the  structural  and  aerodynamic 
properties  of  the  wing. 

2.5  Solving  the  Static  Aeroelastic  System 

Now  that  the  governing  static  aeroelastic  system  of  equations  is  complete,  it  can  be  seen  from 
dimensional  analysis  that  h,  6,  ip2i  and  1^2  are  of  different  dimensions.  Thus,  for  convenience, 
the  system  was  normalized.  A. 2  gives  the  system  properties,  as  well  as  the  quantities  used 
to  normalize  each. 

After  the  system  has  been  normalized,  it  can  be  reduced  to  two  dependent  variables: 
6  and  h.  6  simply  represents  the  twist  about  the  elastic  axis,  which  has  been  assumed  to 
correspond  to  the  reference  axis,  h  represents  the  vertical,  or  plunging,  displacement  of  the 
elastic  axis.  Both  are  now  functions  of  rj,  which  is  the  non-dimensionalized  spanwise  variable, 
X2,  and  runs  from  0  to  1.  This  final  system  can  be  written  as^ 


CHAPTER  2.  PROCEDURE 


25 


-  m4EiQnh"'  +  m^EiQnO"  -  rrii^QnO  +  muQnh'  =  mis  (2.58) 

jy  4mi2  UIqEiQyi  /n  /3  I  /n  A'  z?  /O  u'"  \  Zo  c;n^ 

6* - , - Q  -  mieQnO  +  rrinQnh - BiQ^h  =  mis  +  mig  (2.59) 

m2Ei  +  1  m2Ei  +  1 


With  the  resulting  boundary  conditions: 

?7  =  0: 

miEih'"  —  mim^ElQnh''  +  h'  +  mimsE^O'  =  0  (2.60) 

miEi{m2Ei  +  1)0'”  +  [{m2Ei  -  1)^  +  mim^ElQn]9'  -  mimjElQnh"  =  0  (2.61) 

h  =  e  =  D  (2.62) 

?7  =  1: 

h"  +  msEiQnO  —  m^EiQnh'  =  —m^Ei  (2.63) 

{m2Ei  +  1)9"  +  meEiQnO  -  mrEiQnh'  =  -{ms  +  mg)^!  (2.64) 

h'”  —  miEiQnh”  +  m^EiQnO'  =  0  (2.65) 

{m2Ei  +  1)9'”  —  (4mi2  —  m(,EiQn)9'  —  mlBiQ^h”  =  0  (2.66) 


Here,  Ei  is  the  non-dimensional  transverse  shear  flexibility  for  a  transversely  isotropic 
material,  which  is  defined  as  the  ratio  between  Young’s  modulus  and  the  modulus  of  trans¬ 
verse  shear  rigidity.  If  Kirchoff’s  hypothesis  had  been  assumed,  Ei  would  equal  zero.  The 
m  coefficients  are  functions  of  the  structural  and  aerodynamic  properties  of  the  wing.  A. 3 
gives  expressions  for  the  m  coefficients.  Qn  is  the  normalized  dynamic  pressure.  When 
determining  divergence  speeds,  Qn  will  become  the  eigenvalue  of  the  governing  equations. 
When  deriving  the  wing  mode  shapes,  Qn  will  take  on  a  value  below  the  critical  value  of 
divergence. 
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Figure  2.2:  Simplified  geometric  description  of  wing  plunge,  h,  and  twist,  9. 

Having  reduced  the  system  to  two  dependent  variables,  h  and  6,  which  vary  along  the 
X2  axis.  Fig.  2.1  can  be  simplihed  to  Fig.  2.2.  Notice  that  both  h  and  6  are  descriptions  of 
action  in  the  xi^X:^  plane,  but  they  are  functions  of  X2- 

Due  to  the  fact  that  mode  shapes  become  irrelevant  quantities  after  wing  failure,  the 
hrst  step  in  analysis  is  to  determine  at  what  dynamic  pressure  divergence  occurs.  The 
hrst  step  in  accomplishing  this  is  to  create  a  matrix  inclusive  of  all  wing  properties.  After 
using  MATHEMATICAL  DSolve  function  to  solve  the  governing  equations,  an  8  x  8  matrix. 
Equation  (2.67),  can  be  formed  from  the  boundary  conditions. 
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[A]{/3}  =  {F}  (2.67) 

Where  A  is  the  matrix  describing  the  wing  structural,  geometric,  and  aerodynamic  prop¬ 
erties  dependent  on  airflow,  is  a  vector  made  up  of  the  8  unknown  boundary  conditions, 
and  F  are  the  wing  properties  independent  of  airspeed.  Knowing  that  Qn,  the  normalized 
dynamic  pressure,  is  indicative  of  airspeed  and  density,  it  is  left  as  a  variable.  This  now 
becomes  a  simple  eigenvalue  problem,  and  the  divergence  dynamic  pressure  can  be  obtained 
by  setting  the  determinant  of  the  matrix  equal  to  zero  and  solving  for  Qn. 

In  reality  the  speed  of  MATHEMATICA  forces  the  use  of  a  guess  and  check  technique. 
By  creating  a  loop,  which  inputs  various  values  of  Qn  before  solving  the  governing  equations, 
the  general  trend  can  be  used  to  find  where  the  determinant  equals  zero. 

In  order  to  solve  for  the  bending  shapes,  the  basic  technique  of  matrix  inversion  was  used. 
In  such  a  large  matrix  there  is  the  a  good  chance  of  this  causing  an  ill-conditioned  matrix. 
To  avoid  that  the  process  of  iterative  improvement  was  used. 

{13}  =  |A|-'{F}  (2.68) 


By  using  Equation  (2.68),  the  unknown  boundary  conditions  could  be  determined.  Plug¬ 
ging  these  back  into  the  governing  equations  and  choosing  a  dynamic  pressure  creates  two 
equations,  H{ri)  and  d{ri),  which  describe  the  pitching  and  plunging  motion  of  the  wing  as 
functions  of  the  spanwise  location.  To  generate  truly  useful  data,  these  equations  can  be 
combined  into  a  single  equation  describing  the  effective  angle  of  attack  of  the  wing,  cce//,  as 
follows. 
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Q  _  tan(k)  1^/ 

ae//  =  «o(l  + - (2.69) 

The  MATHEMATIC  A  code  for  both  the  divergence  speed  and  effective  angle  of  attack 
calculations  can  be  found  in  A. 4  and  A. 5,  respectively. 

After  determining  values  for  both  the  divergence  speed  and  effective  angle  of  attack 
distributions  at  different  Mach  numbers,  an  accurate  assessment  of  the  effects  of  transverse 
shear  deformation  in  steady  flow  can  be  made.  The  next  step  in  the  process  is  to  look  at 
unsteady,  oscillatory  flow. 

2.6  Unsteady  Aerodynamics 

While  steady  aerodynamics  follows  basic  principles  of  physics  and  can  be  taught  at  the 
undergraduate  level,  the  topic  of  unsteady  aerodynamics  requires  a  good  deal  more  effort. 
Utilizing  higher  level  math  and  complex  functions,  it  is  an  area  of  study  worthy  of  a  dedicated 
graduate  level  course.  However,  the  unsteady  data  necessary  in  the  context  of  this  analysis 
merely  requires  calculating  values  for  lift  and  moment  from  known  equations.  Somewhat 
complicating  the  process  are  the  Mach  number  effects. 

In  order  to  determine  accurate  values  for  the  aerodynamic  forces,  the  unsteady  equations 
must  be  corrected  for  compressibility.  An  amazingly  complex  process  to  derive,  the  work 
can  be  credited  to  two  men  Theodore  Theodorsen*^^’^®®^  and  P.F.  Jordan^^^’^d  In  the  early 
20th  century,  Theodorsen  derived  a  complex  function,  C{k)  =  F{k)  +iG{k),  used  in  the 
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prediction  of  unsteady  incompressible  aerodynamic  forces,  where  k  is  the  reduced  frequency 
of  oscillation.  This  was  a  giant  leap  for  aerodynamicists  and  mathematicians  at  the  time. 
Shortly  after  the  Second  World  War,  Jordan  took  the  analysis  of  unsteady  flow  a  step  further 
by  correcting  for  the  Mach  number  effects  due  to  compressibility.  Jordan’s  work  directly 
corrects  the  Theodorsen  function  through; 


Ccomp 

-F  i 


comp  ^ 

incomp 

incomp 


where  F 


incomp 


is  the  real  part  of  the  Theodorsen  function  and 


(2.70) 
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In  Equation  (2.71),  the  I  variables  are  known  as  the  Jordan  coefficients  and  can  be  found 
in  Reference  [10].  Once  the  Theodorsen  function  has  been  corrected  for  compressibility,  it 
can  be  used  in  the  equations  for  unsteady  aerodynamics. 

Since  flutter  can  be  described  by  a  neutrally  unstable  oscillation  of  a  wing  at  a  given 
frequency,  u,  lift  and  moment  must  also  be  represented  as  harmonic  functions  with  the  same 
frequency.  That  is: 

L(r/;t)  =  Fe(L(r/)e(*"’*))  (2.72) 

M{r];t)  =  Fe(M(r/)e(*‘"*))  (2.73) 

where  L{ri)  and  M{ri)  are  complex  amplitudes  of  the  unsteady  aerodynamic  loads  and  are 
related  to  the  complex  amplitudes  of  the  oscillatory  modes 


(2,74) 

9{ri;t)  =  Fe(0(r/)e(*"’‘)) 

(2.75) 

by  the  following  relations: 
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L  —  7rpc^u^(hLhh  +  +  hLhh/  +  §Lfi0i)  (2.76) 

M  =  TT  pc^uP'  (JiMoh  +  9Mqq  +  HMoh!  +  dMgo/)  (2.77) 

In  Equations  (2.76-7),  the  coefficients  of  the  mode  shapes  are  functions  of  aerodynamic 
properties  of  the  wing  and  the  Theodorsen  function.  It  is  within  these  coefficients  that  the 
correction  for  compressibility  to  the  Theodorsen  function  is  made.  The  equations  can  be 
found  in  A. 6. 

2.7  Solving  the  Unsteady  Aeroelastic  System 

The  process  for  solving  the  unsteady  aeroelastic  system  is  much  the  same  as  that  used  for  the 
static  state.  Both  are  systems  of  two  governing  equations  and  eight  boundary  conditions. 
However,  a  new  problem  arises  in  the  unsteady  system.  Instead  of  just  one  eigenvalue 
describing  the  static  instability,  there  are  two  for  the  dynamic  instability.  Thus,  instead 
of  dynamic  pressure,  the  flutter  problem  is  a  function  of  hi  and  k,  the  normalized  circular 
frequency  and  reduced  frequency,  respectively.  They  can  be  written  as: 


00 

(2.78) 
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(2.79) 

Aside  from  this  point,  the  unsteady  system  is  similar  to  that  of  the  steady  state.  Sparing 
another  lengthy  derivation,  the  unsteady  model  can  be  written  as: 


+  IT2(s)0(=^)(r/)  =  0 
W3{s)H^^\p)  +  W4{s)e^*\p)  =  0 


(2.80) 

(2.81) 
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with  the  boundary  conditions  at  the  root: 

IT5(s)Lf®(0)  +  IT6(s)0(2)(O)  =  0 

(2.82) 

+  IT8(s)0(=^)(O)  =  0 

(2.83) 

/i(0)  =  0 

(2.84) 

0(0)  =  0 

(2.85) 

and  the  boundary  conditions  at  the  tip; 

Wg{s)H^^\l)  +  Wio{s)e^"\l)  =  0 

(2.86) 

Wu{s)H^^\l)  +  ITi2(s)0®(1)  =  0 

(2.87) 

+  ITi4(s)0«(1)  =  0 

(2.88) 

ITi5(s)if«(l)  +  ITi6(s)0(2)(l)  =  O 

(2.89) 

Here  the  terms  W{s)  is  used  to  describe  the  aerodynamic,  geometric,  and  structural 
properties  of  the  wing,  which  are  not  shown  explicitly.  H  and  0  are  the  Laplace  transforms 
of  h  and  9,  where  the  superscripts  are  used  to  show  the  order  of  each  system. 

Instead  of  using  the  DSolve  command  in  MATHEMATICA,  the  system  was  solved  using 
the  Laplace  transform  technique.  This  involved  taking  the  Laplace  transform  of  the  two 
governing  equations,  factoring  out  the  H  and  0  and  solving  two  equations  for  two  unknowns. 
Once  H  and  0  were  determined  as  functions  of  s,  which  represents  the  independent  variable 
in  the  Laplace  domain,  the  inverse  Laplace  transform  could  be  taken,  giving  h  and  9  as 
functions  of  f].  Now  instead  of  an  eighth  order  system  of  differential  equations,  there  is 
a  system  of  two  algebraic  equations  with  eight  higher  order  unknowns.  This  is  where  the 
Laplace  method  shows  its  true  utility. 
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The  Laplace  Transform  method  takes  into  consideration  a  function’s  initial  conditions 
when  working  with  derivatives.  As  an  example,  the  Laplace  Transform  of  /'  is  sF  —  /(O). 
This  fact  allows  for  the  four  initial  conditions  to  be  taken  into  account  within  the  two 
governing  equations.  Once  this  is  done,  the  system  has  four  unknowns.  By  this  point  in  the 
process  the  two  governing  equations  have  become  quite  large,  and  a  switch  from  pen  and 
paper  derivation  to  a  computerized  solved  is  needed.  Again  MATHEMATICA  hlls  this  role 
quite  nicely. 

Using  MATHEMATICA  and  wing  tip  boundary  conditions,  i.e.  rj  =  1,  a  matrix  can 
be  created  by  substituting  in  the  governing  equations  evaluated  at  the  tip.  This  leaves  a 
four  by  four  matrix  similar  to  the  eight  by  eight  matrix  used  in  the  steady  state  problem. 
The  beneht  of  using  the  Laplace  technique  can  be  seen  when  manipulating  this  matrix,  as 
inverting  it  has  a  signihcantly  lower  chance  of  causing  an  ill-conditioned  matrix. 

This  for  by  four  matrix  now  has  six  unknowns:  the  airflow  properties  k  and  as  well 
as  four  initial  conditions  h'(0),  9'{0),  and  ^"(0).  Depending  on  what  type  of  data  is 

known,  this  system  can  be  used  to  either  determined  flutter  speed  and  frequency  or  the  wing 
mode  shapes. 

In  order  to  determine  flutter  speed  and  frequency,  the  Theodorsen  method  is  used.  This 
method  is  similar  to  the  eigenvalue  problem  used  to  solve  for  divergence  speed,  but  takes 
two  variables  into  account,  k  and  D.  At  hrst  it  would  seem  that  one  matrix  with  two 
variables  could  not  be  solved.  However  since  the  matrix  is  time  dependent,  its  terms  have 
become  complex.  This  means  that  for  the  determinant  to  approach  zero,  both  the  real  and 
imaginary  parts  of  each  term  must  vanish.  Theodorsen  proposed  that  this  could  be  done  by 
selecting  two  values  of  reduced  frequency,  both  close  to  the  assumed  value  that  would  cause 
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the  determinant  to  equal  zero.  Then  the  normalized  frequency  is  varied  until  two  values 
are  determined.  One  where  the  real  part  of  the  determinant  equals  zero,  and  one  where  the 
imaginary  part  equals  zero.  This  process  is  repeated  for  the  other  value  of  k.  The  flutter 
eigenvalues  at  which  the  flutter  determinant  becomes  zero  can  then  be  determined  from  the 
intersection.  This  can  better  be  explained  by  Fig.  2.3. 


Figure  2.3:  Theodorsen’s  method  for  determining  the  flutter  eigenvalues  k  and  O. 

In  this  example,  k  =  0.4  and  0  =  3  would  cause  the  flutter  determinant  to  vanish.  At 
this  point.  Equations  (2.74-5)  can  be  used  to  solve  for  the  flutter  speed,  V,  and  the  flutter 
frequency,  cj. 

Having  now  determined  the  flutter  stability  boundary  corresponding  to  the  maximum 
speed  and  oscillatory  frequency  the  wing  can  sustain,  subcritical  mode  shapes  can  be  deter¬ 
mined.  Utilizing  the  same  method  as  outlined  in  Section  2.5,  the  matrix  can  be  inverted. 
Again  the  process  of  iterative  improvement  reduces  the  chance  of  creating  an  ill-conditioned 
matrix  during  the  inversion  process. 
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2.8  The  Physical  Meaning  of  E/G’ 


Throughout  this  derivation  the  term  E/G’  has  been  included  as  a  structural  property  without 
much  explanation.  However,  it  has  never  been  assigned  real-world  values.  Often  when 
performing  research  it  is  easy  to  focus  on  the  specihcs  and  forget  why  the  research  is  actually 
being  done.  This  situation  can  be  avoided  if  real-life  results  are  always  an  end  goal  of  the 
research. 

In  the  case  of  E/G’,  the  analysis  looks  at  values  from  0-100.  To  better  understand  how 
these  numbers  correlate  to  some  real  world  composits.  Table  5.1  gives  different  materials 
along  with  their  E/G’  values.  Some  materials  such  as  hberglass  and  graphite  are  household 
names.  Some  are  not.  However,  all  of  these  materials  are  used  throughout  industry. 


Material 

E/G’ 

Fiberglass 

7 

Boron 

23 

Graphite 

31 

Table  2.1:  E/G’  Values. 


Chapter  3 

Results  and  Discussion 


Having  laid  out  the  mathematics  behind  the  project,  the  focus  can  now  be  shifted  to  results. 
The  models  developed  in  Chapter  2  are  extremely  generic  models,  which  can  be  used  to 
analyze  a  large  variety  of  wings.  As  long  as  a  wing  hts  into  a  fairly  general  category, 
its  specihc  properties  can  be  input  into  the  model  and  results  for  critical  eigenvalues  and 
subcritical  mode  shapes  can  be  determined.  For  this  analysis  the  wings  described  in  above 
will  be  used  to  analyze  the  effects  of  transverse  shear  rigidity  on  the  static  and  dynamic 
aeroelastic  instabilities  in  the  presence  of  compressibility  effects. 

3.1  Wing  Models 

Once  a  generic  model  was  developed,  it  became  necessary  to  gather  information  on  wings 
of  known  geometric,  aerodynamic,  and  structural  properties.  This  is  necessary  in  order  to 
have  specihc  wings  to  analyze. 
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Wing 

Goland’s 

400R 

Span  (ft) 

40.0 

2.90 

Ghord  (ft) 

6.0 

0.79 

AR 

6.67 

3.69 

Thickness 

0.1 

0.04 

Table  3.1:  Goland’s  Wing  and  the  400R  Wing. 

The  two  wings  chosen  were  Goland’s  wing^^  and  the  400R^^  wing.  Specihc  properties  of 
both  can  be  found  in  Table  3.1.  Obviously,  the  400R  wing  is  much  smaller  and  was  most 
likely  designed  for  wind  tunnel  testing.  Looking  at  the  non-dimensional  properties  of  aspect 
ratio  and  thickness,  it  can  be  seen  that  the  400R  wing  is  much  shorter  and  thinner  than 
Goland’s  wing.  In  terms  of  this  analysis,  as  a  structure  become  thicker  it  tends  to  be  more 
affected  by  transverse  shear  deformation. 

Unfortunately  neither  of  these  wings  has  a  very  low  transverse  shear  moduli  causing 
E/G'  ~  0.  Thus  both  are  considered  to  be  made  of  classical  materials.  This  problem 
is  overcome  by  simply  varying  the  value  of  E/G'.  As  long  as  the  basic  geometric  and 
aerodynamic  properties  of  the  wings  are  held  constant,  varying  only  a  few  properties  allows 
for  a  parametric  analysis  to  be  performed  on  the  aeroelastic  instability  of  the  wing. 

3.2  Test  Equipment 
3.2.1  MATHEMATICAL 

All  computer  work  was  done  using  the  Silicon  Graphics  system  in  the  USNA  GADIG  center. 
A  UNIX  based  version  of  MATHEMATIC  A  was  run.  Input  could  be  made  via  a  TELNET 
connection  from  outside  of  the  facility.  However,  direct  input  into  one  of  the  GADIG  servers 
was  more  efficient. 
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Although  a  powerful  computational  tool  and  the  choice  of  many  mathematicians  for 
symbolic  manipulation,  MATHEMATIC  A  was  shown  to  have  some  flaws.  The  main  problem 
encountered  was  the  black  box  design  of  the  program.  Even  though  MATHEMATIC  A  has 
literally  thousands  of  functions  available  for  symbolic  manipulation  and  numeric  analysis, 
very  few  allow  the  user  insight  into  how  the  functions  actually  work.  While  this  is  acceptable 
for  basic  math  functions,  processes  such  as  matrix  inversion  and  the  numerical  solver  may 
or  may  not  be  operating  as  the  user  desires. 

3.3  Steady  Flow  Analysis 

3.3.1  Effect  of  Mach  Number  on  Divergence  Speed 

Before  any  meaningful  analysis  can  be  performed  on  the  what  shape  a  wing  will  take  at  a  give 
air  speed,  it  is  necessary  to  determine  at  what  speed  the  wing  will  fail.  This  is  accomplished 
by  the  method  laid  out  in  Section  2.5  with  the  Goland  wing  properties.  For  this  analysis 
the  wing  will  be  swept  forward  20  deg.  By  varying  Mach  number  and  the  transverse  shear 
rigidity,  a  set  of  data  points  was  created,  which  showed  the  effects  of  both  on  the  divergence 
speed.  Fig.  3.1  shows  these  results. 

Divergence  spped  is  measured  as  the  nondimensional  value,  Qn,  which  is  a  function  of 
airspeed  and  air  density.  Notice  the  decrease  in  divergence  speed  as  Mach  number  increases. 
This  is  representative  of  the  higher  amount  of  kinetic  energy  carried  in  high  speed  flow.  It  is 
important  to  remember  that  Mach  number  can  be  used  to  measure  how  compressed  the  air 
has  become.  When  the  freestream  air  is  compressed  by  the  wing  at  higher  Mach  numbers,  it 
will  transfer  a  greater  amount  of  energy  into  the  structure.  The  wing  on  the  other  hand  can 
only  absorb  a  constant  amount  of  this  in  the  form  of  strain  energy.  Thus,  for  higher  Mach 
number  divergence  occurs  at  a  lower  value  of  Qn- 
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Figure  3.1:  Critical  dynamic  pressure  vs.  Mach  number  for  various  values  of  transverse  shear 
rigidity. 

One  may  wonder  what  signihcance  Fig.  3.1  has  knowing  that  both  dynamic  pressure  and 
Mach  number  are  related  to  airspeed.  To  add  atmospheric  relevance,  the  two  lines  drawn 
in  gray  represent  the  range  of  values  possible  in  the  Earth’s  atmosphere.  The  upper  line 
represents  sea  level  conditions  on  a  standard  day.  The  lower  line  represents  50,000  feet  on  a 
standard  day.  This  shows  that  at  sea  level  higher  dynamic  pressures  can  be  absorbed  by  the 
wing,  while  at  50,000  feet  higher  Mach  numbers  can  be  attained  before  divergence  occurs. 
These  results  represent  another  beneht  gained  by  high  performance  aircraft  at  high  altitudes. 
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3.3.2  Effects  of  Compressibility  on  Effective  Angle  of  Attack 

Having  determined  the  divergence  speeds  for  the  Goland  wing  at  various  atmospheric  con¬ 
ditions  and  Mach  numbers,  analysis  of  the  wing  deformation  can  begin.  Again  using  the 
Goland  wing,  the  effective  angle  of  attack  along  the  wing  semi-span  was  determined  at  var¬ 
ious  airspeeds  relating  approximately  to  Mach  0.3,  0.7,  and  0.8  at  sea  level.  These  results 
are  shown  in  Fig.  3.2. 


Figure  3.2:  Spanwise  distribution  of  effective  angle  of  attack  across  wing  semi-span. 

The  first  point  to  be  made  concerns  the  effects  of  compressibility  at  low  airspeeds.  Fig. 
3.2a  shows  that  at  200  knots,  the  limit  of  incompressible  flow  theory,  the  compressible 
analysis  actually  shows  a  smaller  effective  angle  of  attack  distribution  across  the  wing.  Even 
though  the  difference  is  almost  negligible,  a  fraction  of  a  degree  at  the  tip,  this  adds  a  factor 
of  safety  to  the  design  of  low  speed  aircraft  whose  designers  use  incompressible  flow  theory. 


CHAPTER  3.  RESULTS  AND  DISCUSSION 


40 


For  instance,  a  general  aviation  aircraft  intended  for  200  knots,  approximately  Mach  0.3  at 
sea  level,  was  most  likely  designed  using  incompressible  flow  theory.  When  the  compressible 
analysis  is  performed  in  this  performance  region,  it  shows  that  there  is  actually  less  deforma¬ 
tion  than  previously  calculated.  As  general  aviation  aircraft  companies,  such  as  Cessna  and 
Diamond,  move  to  more  and  more  composite  aircraft  components,  this  information  could 
become  very  useful. 

Another  interesting  occurrence  is  the  divergence  speed.  Fig.  3.1  showed  that  for  the 
Goland  wing  divergence  should  occur  between  Mach  0.7  and  0.8.  Fig.  3.2  shows  this  hap¬ 
pening.  While  the  jump  between  an  airspeed  of  200,  Fig.  3.2a,  and  460  knots.  Fig.  3.2b, 
only  caused  an  8  degree  change  in  the  effective  angle  of  attack  at  the  wing  tip,  the  jump 
between  460  and  530  knots.  Fig.  3.2c,  caused  a  larger  jump  in  both  the  compressible  and 
incompressible  analysis.  This  reinforces  the  fact  that  for  the  Goland  wing,  divergence  speed 
does  occur  in  this  range. 

Forward  swept  wings  tend  to  have  a  lower  divergence  speed  than  straight  or  swept  back 
wings.  Forward  swept  wings  also  show  more  deformation  than  straight  or  swept  back  wings 
at  the  same  airspeed.  This  fact  was  a  major  problem  when  designing  the  X-29  Forward 
Swept  Wing  aircraft.  To  alleviate  this  problem,  composite  tailoring  was  used.  In  contrast, 
swept  back  wings  experience  divergence  at  much  higher  speeds  than  swept  forward  wings, 
therefore  lessening  wing  deformation.  Fig.  3.3  shows  the  effective  angle  of  attack  across  the 
wing  semi-span  for  Goland’s  wing  swept  back  20  deg.  Results  are  shown  for  airspeeds  of  460 
and  530  knots,  relating  to  Mach  0.7  and  0.8  at  sea  level. 
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Figure  3.3:  Effective  angle  of  attack  across  semi-span  for  Goland’s  wing  swept  back  20deg. 

Using  Fig.  3.3  it  is  easier  to  point  out  how  little  changes  in  airspeed  affect  wing  defor¬ 
mation.  While  there  is  a  20  deg  to  30  deg  jump  between  460  and  530  knots  for  the  forward 
swept  wing,  there  is  only  a  0.5  deg  jump  for  the  swept  back  wing. 

3.3.3  Effect  of  Sweep  on  Effective  Angle  of  Attack 

Although  it  is  quite  simple  to  state  that  forward  swept  wings  are  deformed  more  than  swept 
back  wings  at  similar  air  speeds,  further  investigation  into  the  effects  of  sweep  produce 
interesting  results.  Before  these  results  can  be  discussed  it  is  necessary  to  explain  the  level 
of  influence  sweep  has  on  the  aerodynamic  and  geometric  properties  of  a  wing. 
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When  air  flows  across  a  wing,  only  the  normal  component  of  the  airspeed  to  the  leading 
edge  of  the  wing  creates  aerodynamic  forces  and  moments.  Thus,  if  an  aircraft  with  wings 
swept  at  an  angle  A,  is  flying  at  an  airspeed  V,  then  the  effective  airspeed  is  actually  only 
V  cos(A).  This  relationship  can  be  directly  applied  to  Mach  number.  In  truth  if  an  aircraft 
is  flying  at  Mach  1  with  its  wings  swept  back  20  deg,  its  wing  cross-section  is  only  seeing 
Mach  0.94.  Sweep  angle  also  has  an  effect  on  the  three  dimensional  lift-curve  slope  of  a 
wing,  oq.  As  the  magnitude  of  sweep  increases,  oq  also  increases  by  1/ cos(A).  This  effect 
acts  opposite  to  the  loss  of  lift  caused  by  sweep’s  reduction  of  effective  airspeed. 

Geometrically,  sweep  alters  the  effective  angle  of  attack  as  shown  in  Equation  (2.56). 
tte//  is  a  coupling  of  the  bending  and  torsional  deformations,  otherwise  known  as  plunging 
and  pitching  modes.  Sweep  angle  is  necessary  to  relate  these  two.  However,  instead  of 
using  the  cosine  of  sweep,  this  coupling  relies  on  the  tangent  of  sweep.  Knowing  that  the 
tangent  function  diverges  at  90deg  and  -90deg,  special  attention  must  be  paid  to  its  effect  on 
the  results.  Obviously  as  sweep  approaches  these  values,  the  two-dimensional  aerodynamic 
theory  used  in  the  model  will  no  longer  work. 

With  all  these  sweep  effects  working  with  and  against  each  other,  as  both  cosine  and 
tangent  functions,  it  is  hard  to  predict  exactly  what  overal  effect  sweep  will  have.  Fortunately, 
computers  can  perform  these  complicated  calculations  in  seconds  and  output  the  results. 
From  that  Fig.  3.4  shows  the  effects  of  sweep  on  effective  angle  of  attack.  Positive  values  of 
A  indicate  a  swept  back  wing  and  negative  values  indicate  a  forward  swept  wing. 
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Figure  3.4:  Effect  of  sweep  on  effective  angle  of  attack  at  550  knots. 


At  first  glance,  the  results  for  the  swept  back  appear  reversed.  Not  only  is  the  magnitude 
smaller,  but  the  wing  actually  has  a  lower  angle  of  attack  at  the  tip  than  the  root.  This 
is  due  to  tan(A)  being  positive  for  a  swept  back  wing.  Even  though  the  twist  pulls  the 
wing  tip  up,  the  coupled  effect  of  the  wing  plunging  mode  causes  an  overall  decrease  in  the 
effective  angle  of  attack.  Eor  a  swept  back  wing  these  two  offset  one  another,  but  in  swept 
forward  wings  they  both  work  together.  This  likely  explains  the  fact  that  forward  swept 
wings  diverge  at  lower  air  speeds. 
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The  truly  interesting  results  come  from  the  different  values  of  forward  sweep.  Start¬ 
ing  from  10  deg,  increasing  forward  sweep  causes  an  increase  in  effective  angle  of  attack. 
However,  this  trend  reverses  somewhere  between  A  =  SOdeg  and  A  =  50deg.  By  the  time 
A  =  60deg  the  deformation  is  back  below  A  =  2Ddeg.  This  occurrence  is  most  likely  due  to 
the  switch  between  a  cosine  dominant  function  and  a  tangent  dominated  function.  At  -30 
deg  cosine  is  0.87  and  tangent  is  -0.58,  by  -50  deg  cosine  is  0.64  while  tangent  is  1.19. 

This  switch  between  dominant  terms  in  Equation  (2.56)  signihes  two  points.  The  hrst 
is  that  at  a  certain  sweep  angle,  geometric  properties  become  more  important  than  purely 
aerodynamic  principles.  The  second,  which  is  less  encouraging,  is  that  at  a  certain  point,  the 
equations  may  become  unrealistic.  To  validate  results  at  this  transition  point  experimental 
data  would  have  to  be  collected.  If  indeed  this  data  did  not  show  this  region  of  transition, 
then  a  line  would  have  to  be  drawn  showing  where  the  theory  no  longer  held  true. 

3.3.4  Effects  of  Transverse  Shear  Deformation  on  Effective  Angle 
of  Attack 

Having  analyzed  the  effects  of  compressibility  and  sweep  on  the  steady  state  aeroelastic 
properties  on  Goland’s  wing,  an  accnrate  discussion  of  the  effects  of  transverse  shear  can 
now  be  made.  It  was  already  shown  in  Section  5.1.1  that  as  E/G’  increased,  divergence  speed 
decreased.  This  made  sense  as  more  flexibility  throngh  the  wing’s  thickness  would  accelerate 
wing  instability.  Now  comparing  the  effect  of  transverse  shear  deformation  on  effective  angle 
of  attack,  more  in  depth  conclusions  can  be  made.  Fig.  3.5  shows  the  spanwise  effective  angle 
of  attack  distribution  for  various  values  of  E/G’  for  both  incompressible  and  compressible 
flows. 


CHAPTER  3.  RESULTS  AND  DISCUSSION 


45 


Figure  3.5:  Effect  of  transverse  shear  on  effective  angle  of  attack  at  200  knots. 

Again  similar  to  the  results  drawn  from  the  analysis  of  divergence  speeds,  it  can  be  seen 
from  the  analysis  of  Fig.  3.5  that  the  more  flexible  the  wing  in  transverse  shear,  the  more  it 
will  deform  at  a  certain  air  speed.  This  can  directly  be  related  to  a  lower  divergence  speed 
for  flexible  wings. 

Also  from  Fig.  3.5a  it  can  be  seen  that  at  200  knots,  incompressible  analysis  is  fairly 
similar  for  all  valnes  of  E/G’.  There  is  less  than  a  1  %  difference  in  the  effective  angle  of 
attack  at  the  wing  tip.  However,  the  compressible  analysis  shows  that  a  mnch  larger  gap 
in  bending  occurs  between  the  different  valnes  of  E/G’.  An  exact  physical  explanation  for 
this  is  difhcnlt  to  draw.  One  explanation  may  be  that  as  the  wing  bends,  compressible  flow 
has  a  high  tendency  to  increase  the  aerodynamic  loads  and  moments  it  places  on  the  wing. 
Thus,  a  change  in  wing  flexibility  would  have  a  greater  effect  in  compressible  flow,  than  in 
incompressible  flow.  It  is  important  to  remember  that  these  resnlts  and  hypothesis  are  for 
low  speed  flows.  Because  as  Fig.  3.6  points  ont,  varying  three  properties  makes  an  analysis 
of  results  very  difficult. 
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Figure  3.6:  Effect  of  transverse  shear  on  effective  angle  of  attack  at  460  knots. 


Fig.  3.6  shows  that  even  at  high  airspeeds  incompressible  theory  causes  a  larger  deforma¬ 
tion  than  compressible  theory.  However,  the  gap  is  dehnitely  closing.  At  valnes  of  E/G'  =  0 
there  is  hardly  any  difference  between  the  compressible  and  incompressible  analysis.  It  is 
qnite  possible  that  at  higher  airspeeds  and  thus,  higher  Mach  numbers,  the  magnitude  of  the 
compressible  aerodynamic  loads  begins  to  catch  np  with  that  of  the  incompressible  loads. 

To  further  back  up  this  argument,  the  large  deformation  gap  that  was  seen  in  the  com¬ 
pressible  analysis  at  200  knots  has  been  signihcantly  rednced.  Fig.  3.6  actually  shows 
that  the  incompressible  theory  now  has  a  larger  gap  between  the  wing  tip  deformations  of 
E/G'  =  0  and  E/G'  =  100.  To  explain  this  phenomenon  better,  deeper  research  into  com¬ 
pressible  snbsonic  aerodynamics  wonld  be  necessary.  As  these  results  seem  to  go  against 
intnition,  it  appears  that  this  is  simply  a  case  of  Ending  answers,  which  lead  to  more  qnes- 
tions. 
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3.4  Flutter  Analysis 

3.4.1  The  Influence  of  Transverse  Shear  Deformation  on  Flutter 
Frequency 

Transitioning  now  to  the  unsteady  side  of  the  held,  the  focus  switches  from  divergence 
speeds  and  effective  angles  of  attack  to  hutter  frequencies  and  velocities.  It  is  important  to 
remember  throughout  this  analysis  that  wing  failures  can  occur  through  both  static  aeroe- 
lastic  instability,  divergence,  and  dynamic  instability,  hutter.  In  diherent  conhgurations  and 
environments,  either  type  of  failure  could  occur  hrst.  Thus,  determining  the  conditions  nec¬ 
essary  for  both  categories  to  occur  is  vital  in  understanding  what  how  regimes  an  aircraft 
can  operate  in. 

Moving  on,  Theodorsen’s  method  was  utilized  to  investigate  the  ehects  of  transverse  shear 
rigidity  and  Mach  number  on  the  hutter  eigenvalues.  Unlike  the  static  instability  described 
by  a  single  eigenvalue,  the  hutter  instability  is  described  by  two  eigenvalues.  Both  of  these 
eigenvalues  were  normalized  during  this  analysis  in  order  to  more  eliminate  the  necessity  of 
carrying  units.  As  hutter  frequency,  rate  of  wing  oscillation  at  failure,  is  seemingly  more 
difficult  to  understand  it  will  be  dealt  with  hrst.  Fig.  3.7  shows  the  ehects  of  transverse 
shear  rigidity. 
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Figure  3.7:  Effect  of  transverse  shear  flexibility  on  flntter  freqnency. 


From  this  it  can  be  seen  that  as  the  wing  becomes  more  flexible,  its  flntter  frequency 
decreases  at  all  Mach  numbers.  However,  as  the  Mach  number  increases  up  into  the  high 
speed  subsonic  and  snpersonic  regions  the  curves  become  more  linear.  This  follows  the 
hypothesis  from  steady  flow  that  transverse  shear  rigidity  has  a  larger  impact  npon  wings  in 
low  speed  snbsonic  flow.  It  is  also  clear  that  as  the  wing  becomes  more  flexible  in  transverse 
shear  flutter  occurs  at  lower  frequencies.  This  phenomenon  is  consistant  with  the  inability 
of  the  wing  structnre  to  dissipate  the  energy  absobed  due  to  the  aerodynamic  loads. 
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Figure  3.8:  Effect  of  Mach  number  on  flutter  frequency. 

To  better  illustrate  this  idea,  Fig.  3.8  displays  flutter  frequency  against  Mach  number 
for  three  values  of  E/G’.  If  Fig.  3.7  left  any  doubt  that  transverse  shear  deformation  is  a 
problem  at  low  airspeeds,  Fig.  3.8  will  erase  it. 

The  results  show  the  importance  of  Mach  number  on  flutter  frequency.  Notice  the  large 
difference  between  flutter  frequencies  at  Mach  0.  The  wing  with  E/G'  =  10  has  a  flutter 
frequency  approximately  20%  larger  than  the  E/G'  =  100  wing.  As  the  Mach  number 
increases,  this  difference  decreases.  At  Mach  1  this  difference  has  fallen  to  only  12%,  and 
at  Mach  2  it  is  well  under  10%.  Again,  the  effects  of  transverse  shear  deformation  on  the 
flutter  frequency  are  greater  at  the  lower  Mach  numbers. 
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Fig.  3.8  also  draws  an  excellent  picture  of  what  is  commonly  known  as  ’’transonic  dip”. 
Notice  the  sudden  decrease  in  the  flutter  frequency  as  the  Mach  number  crosses  from  0.95 
to  1.2.  Past  Mach  1.2  the  curve  again  levels  out.  This  occurrence  is  common  in  transonic 
aerodynamic  and  gives  shows  theoretically  just  how  dynamic  the  transonic  regime  can  be. 

3.4.2  The  Influence  of  Transverse  Shear  Deformation  on  Flutter 
Speed 

Flutter  speed,  the  velocity  at  which  dynamic  instability  occurs,  can  be  dealt  with  much  like 
flutter  frequency.  Exhibiting  many  of  the  same  tendencies  when  transverse  shear  rigidity 
and  Mach  number  are  altered,  the  flutter  speed  is  much  simpler  to  physically  comprehend. 
Basically  above  a  certain  airspeed,  an  oscillating  wing  will  no  longer  be  able  to  dissipate  the 
aerodynamic  energy  absorbed  by  the  wing  structure.  When  the  structure  absorbs  energy  at  a 
rate  higher  than  it  can  dissipate  that  energy,  wing  oscillations  will  increase  in  amplitude  until 
structural  failure  occurs.  The  flutter  speed  is  the  neutrally  stable  speed,  which  separates 
damped  oscillations  from  undamped  oscillations. 

In  the  test  flights  of  many  aircraft  certain  ’’danger  zones”  are  encountered  when  the 
flutter  speed  is  reached.  However,  often  higher  air  speeds  can  be  attained  simply  by  changing 
altitude  and  circumventing  the  danger  zone.  If  a  change  in  altitude  can  be  used  to  avoid 
certain  flutter  speeds  it  is  not  unlikely  that  Mach  effects  have  something  to  do  with  flutter 
speeds.  Fig.  3.9  shows  the  effects  of  tranverse  shear  flexibility  on  flutter  speed  for  various 
Mach  numbers. 
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Figure  3.9:  Effect  of  transverse  shear  on  flutter  speed. 

Similar  to  its  effect  on  flutter  frequency,  an  increase  in  transverse  shear  tends  to  lower 
the  flntter  speed.  It  is  interesting  to  note,  that  as  E/G’  becomes  large,  the  Mach  lines  seem 
to  converge  upon  one  another.  All  subsonic  lines  seem  to  cnrve  np  towards  an  imaginary 
tangential  line,  while  the  snpersonic  line  curves  down  towards  the  same  line.  It  turns  out  that 
this  line  seems  to  be  exactly  where  the  Mach  1  results  would  be  depicted.  Unfortunately, 
aerodynamic  data  for  Mach  1  is  extremely  unstable  and  thus  this  point  cannot  be  dehnitively 
proven.  However,  drawing  a  straight  line  from  where  Mach  1  should  cross  the  dependent  axis 
to  the  convergence  of  the  other  four  curves,  it  does  not  seem  unreasonable  that  the  curves 
could  be  converging  to  Mach  1. 

In  gas  dynamics  both  subsonic  and  snpersonic  flows  tend  to  approach  Mach  1  at  the 
throat  of  nozzles.  Could  it  be  that  flutter  speed  exhibits  the  same  properties?  It  looks  as 
if  both  the  subsonic  and  supersonic  curves  are  converging  towards  the  imaginary  Mach  1 
line  as  E/G’  increases.  Fnrther  investigation  of  this  behavior  could  produce  very  interesting 
results. 
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Getting  back  to  the  effects  of  Mach  number  on  flutter  speed,  Fig.  3.10  shows  once  again 
the  inverse  effect  of  Mach  number  on  aeroelastic  properties.  As  the  Mach  number  increases, 
the  difference  between  flutter  speeds  for  various  values  of  E/ G’  decreases.  Similar  to  Fig.  3.1, 
the  grey  atmospheric  boundary  lines  have  been  placed  on  the  plot.  The  left  line  represents 
sea  level  conditions,  while  the  right  line  represents  50,000  ft.  Notice  that  for  flutter  speed, 
the  atmospheric  range  is  much  smaller  than  for  divergence  speed. 


Figure  3.10:  Effect  of  Mach  number  on  flutter  speed. 


As  mentioned  before,  the  three  curves  converge  as  Mach  number  is  increased  up  to  a 
point.  Fig.  3.10  actually  shows  a  slight  divergence  between  the  curves  past  Mach  1.2.  This 
could  be  due  to  the  E/G'  =  10  curve  showing  somewhat  erratic  behavior,  or  possibly  a  new 
phenomenon  not  yet  seen. 
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The  inclusion  of  the  atmospheric  boundary  lines  shows  where  possible  inaccuracies  could 
have  occurred  in  a  merely  incompressible  analysis.  By  the  time  the  flutter  speed  curves  reach 
the  hrst  boundary  line,  sea  level,  they  have  already  decreased  almost  10%.  Similar  to  the 
trend  shown  in  the  divergence  speed  analysis,  the  flutter  speed  decreases  as  Mach  number 
increases.  This  could  be  due  to  the  increased  energy  in  the  compressed  flow. 

3.4.3  Comparison  Between  Goland’s  Wing  and  the  400R  Wing 

Having  looked  at  what  effects  both  Mach  number  and  transverse  shear  flexibility  have  on 
the  aeroelastic  properties  of  Goland’s  wing,  the  next  step  is  to  look  into  the  effects  that  wing 
geometry  has  on  the  flutter  eigenvalues.  As  shown  in  Table  3.1,  the  400R  wing  is  a  thinner 
wing  than  Goland’s  wing.  It  also  has  a  lower  aspect  ratio.  Fig.  3.11  shows  the  comparison 
between  the  flutter  frequencies  for  both  Goland’s  wing  and  the  400R  wing. 

Notice  that  as  the  tranverse  shear  flexibility  is  increased,  the  flutter  frequency  of  Goland’s 
wing  is  more  affected.  Stepping  back  a  moment  to  think  about  this,  the  answer  to  why 
this  occurs  becomes  apparent.  Transverse  shear  is  a  measure  of  the  deformation  occurring 
through  the  thickness  of  the  wing.  If  one  wing  is  thicker  than  another,  such  as  is  the  case 
with  Goland’s  wing  and  the  400R,  it  would  make  sense  that  the  thicker  wing  is  more  affected 
by  a  change  in  transverse  shear  flexibility. 
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Figure  3.11:  Comparison  of  flutter  frequencies  for  Goland’s  wing  and  the  400R. 

The  effects  of  wing  thickness  on  flutter  speed  produce  much  the  same  results.  Fig.  3.12 
shows  the  flutter  speed  comparison. 

Again  an  increase  in  the  value  for  E/G’  produces  a  noticeable  change  in  the  flutter 
speed  for  Goland’s  wing.  However,  almost  no  change  occurs  to  the  characteristics  of  the 
400R.  Again  this  should  concern  general  aviation  enthusiasts  more  than  high  performance 
designers.  Most  low  speed  aircraft  utilize  a  thick  wing  to  generate  lift  at  low  airspeeds.  High 
performance  jets  can  create  most  of  their  lift  through  airspeed  and  require  small  thickness 
and  camber  in  their  wings.  Again  it  turns  out  that  transverse  shear  deformation  affects 
general  aviation  aircraft  more  than  their  high  performance  counterparts. 
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Figure  3.12:  Comparison  of  flutter  speeds  for  Goland’s  wing  and  the  400R. 

3.5  Summary  and  Conclusions 

Initially,  this  research  was  focused  on  analyzing  the  aeroelastic  effects  of  compressibility  and 
transverse  shear  deformation  with  the  thought  that  the  hndings  could  be  used  to  design  high- 
performance  aircraft.  However,  as  the  results  have  shown  it  is  quite  obvious  that  transverse 
shear  is  a  larger  problem  at  low  air  speeds.  Every  analysis  from  the  effect  of  transverse  shear 
flexibility  to  wing  thickness  showed  that  the  wing  properties  common  in  general  aviation 
aircraft  led  to  larger  variations  in  aerodynamic  performance. 
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This  is  not  to  say  that  high  performance  aircraft  should  have  no  concern  of  aeroelastic 
failure.  It  was  shown  that  as  Mach  number  increases  all  critical  speeds  decrease,  which  works 
against  the  jet  community.  Nevertheless,  as  all  aviation  communities  from  ’’fighter  jocks”  to 
civilian  student  pilots  continue  the  search  for  lighter,  faster,  and  more  dependable  aircraft, 
composite  structures  will  no  doubt  be  used  in  place  of  their  metallic  ancestors. 

It  is  important  to  remember  that  although  theoretical  analysis  has  its  flaws,  it  is  useful  as 
a  low-cost  method  of  determining  trends  in  data.  Specific  trends  discerned  from  this  analysis 
are: 

a.  At  low  air  speeds,  and  thus  low  Mach  numbers,  variations  in  transverse  shear  flexibility 
had  more  of  an  effect  on  a  wing’s  aeroelastic  instabilities  than  at  high  speeds. 

b.  As  Mach  number  increases,  the  critical  air  speeds  in  both  steady  and  unsteady  flow 
decrease. 

c.  Compressible  analysis  varies  from  incompressible  analysis  by  anywhere  from  0-30%  in 
standard  atmospheric  conditions. 
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Table  A.l:  Definitions  of  aerodynamic  characteristics. 
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Table  A. 2:  Non-dimensionalized  values. 


V  =  f 

d  1  d 

dx2  I  drj 

p  _  1 

dx2  drj'^ 


e  = 


e  =  e 


Xo  =  ^ 


Ei  = 


G' 


ARi  =  ^ 


Qn 


ARl(l-u^) 
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Table  A. 3:  Determination  of  ‘m’  coefficients. 

"^2  =  ^  rni2  =  ^ 

7713  =  mittoARi  mi3  =  ^ 

=  rriiaotanA  mi4  =  ^ 

ms  =  miOottoAi^iQ^  mis  =  ^ 

"^6  =  "^16  =  ^ 

™  ™  tanA  ^  my 

"^7  -  me^  mi7  -  — 

ms  =  meaoQn  mis  =  ^ 
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Table  A. 4:  MATHEMATICA  code  for  Divergence  Speed. 


Qnlist  =  {}; 

%  Dehne  nullset  matrix 

Do[ 

%  Begin  loop  to  cycle  through  values  of  Qm 

pois  =  .25; 

%  Set  Poison’s  ratio 

El  =  70; 

%  Set  E/G’ 

Clear  [h]; 

Clear  [theta]; 

Clear  [equations]; 

%  Clear  all  variables 

Clear  [A]; 

Clear  [sol]; 

sweep  =  -.3509; 

%  Sweep 

aO  =  6. 28/(1  -|-  4/8  Cos[sweep]); 

%  Incompressible  lift  curve  slope 

AR  =  8; 

%  Aspect  ratio 

ARl  =  AR/2; 

%  Semi-span  aspect  ratio 

alphaO  =  0.0873; 

%  Angle  of  attack 

alphabar  =  alphaO; 

%  Angle  of  attack  (1st  normalization) 

alphadbar  =  alphaO; 

%  Angle  of  attack  (2nd  normalization) 

ebar  =  .1; 

%  Elastic  offset 

tbar  =  .1; 

%  Wing  thickness 

cmac  =  0; 

%  Moment  coefficient 

ml  =  (tbar)^/(12  (1  -  (pois)^)  (ARl)^); 

m2  =  (tbar) ^/ (2  (1  -|-  pois)); 

m3  =  ml  aO  ARl; 

m4  =  ml  aO  Tan  [sweep]; 

m5  =  ml  aO  alphabar  ARl  Qn; 

%  Dehne  coefficients  from  governing  equations 

m6  =  ((tbar)^/(ARl  (1  -  (pois)^)))*ebar  aO; 

m7  =  m6  Tan  [sweep] /ARl; 

m8  =  m6  alphadbar  Qn; 

m9  =  ((tbar)^/((l  -  (pois)^)  ARl))*cmac  Qn; 
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ml2  =  m2/ml; 

ml3  =  m3/ml; 

ml4  =  m4/ml; 

ml5  =  m5/ml; 

%  Define  coefficients  from  governing  equations 

ml6  =  m6/ml; 

ml7  =  m7/ml; 

ml8  =  m8/ml; 

ml9  =  m9/ml; 

al  =  -m4  El  Qn; 

a2  =  m3  El  Qn; 

a3  =  -ml3  Qn; 

%  Define  coefficients  from  governing  equations 

a4  =  ml4  Qn; 

bl  =  ml5; 

a5  =  -(4  ml2  -  m6  El  Qn)/(m2  El  +  1); 
a6  =  -ml6  Qn; 

a7  =  ml 7  Qn;  %  Define  coefficients  from  governing  equations 

a8  =  -(m7/(m2  El  +  1))  El  Qn; 
b2  =  ml8  +  ml9; 

cl  =  ml  El; 

c2  =  -ml  m4  El^  Qn; 

c3  =  1; 

c4  =  ml  m3  El^  Qn;  %  Define  coefficients  from  governing  equations 

c5  =  ml  El  (m2  El  -|-  1); 
c6  =  ((m2  El  -|-  1)^  -|-  ml  m6  El^  Qn); 
c7  =  -ml  m7  El^  Qn; 

dl  =  m3  El  Qn; 
d2  =  -m4  El  Qn; 
b3  =  -m5  El; 
d3  =  (m2  El  +  1);; 
d4  =  m6  El  Qn; 
d5  =  -m7  El  Qn; 
b4  =  -(m8  -|-  m9)  El; 
d6  =  -m4  El  Qn; 
d7  =  m3  El  Qn; 
d8  =  (m2  El  +  1); 
d9  =  -(4  ml2  -  m6  El  Qn); 
dlO  =  -m7  El  Qn; 


%  Define  coefficients  from  governing  equations 
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sol  =  DSolve[ 

{  h””[y]  +  al  h”’[y]  +  a2  theta”  [y] 

+  a3  theta[y]  +  a4  h’[y]  ==  bl, 
theta”  ”  [y]  +  a5  theta”  [y] 

+  a6  theta[y]  +  a7  h’[y]  +  a8  h”’[y]  == 
b2},  {h[y],  theta[y]},  y]; 


h[y_]  =  First [h[y]  /.  sol]; 
theta[y_]  =  First  [theta  [y]  /.  sol]; 
equations  =  {h[0]  ==  0, 
theta[0]  ==  0, 
cl  h”’[0]  +  c2  h”[0] 

+  c3  h’[0]  +  c4  theta’ [0]  ==  0, 
c5  theta” ’[0]  +  c6  theta’ [0] 

+  c7  h”[0]  ==  0, 

h”[l]  +  dl  theta[l]  +  d2  h’[l]  ==  b3, 
d3  theta”  [1]  +  d4  theta[l]  , 

+  d5  h’[l]  ==  b4 

h’”[l]  +  d6  h”[l]  +  d7  theta’ [1]  ==  0, 
d8  theta” ’[1]  +  d9  theta’ [1] 

+  dlO  h”[l]  ==  0}  ; 


equations  =  Simplify  [equations]; 

A  =  Table  [Coefficient  [equations  [[i,  1]], 
C[j]],  {i,  1,  8},  {j,  1,  8}]; 


detA  =  Det[A]; 

Print[detA]; 

Qnlist  =  Append [Qnlist,  detA], 
{Qn,  X,  Y,  Z}] 


%  Solve  governing  equations 


%  Input  boundary  conditions 


%  Simplify  system 
%  Create  matrix  of  coefficients 


%  Calculate  determinant 
%  Output  determinant 
%  Create  list  of  determinants 
%  Run  loop  from  X  to  Y  by  Z 
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Table  A.5:  MATHEMATICA  code  for  Effective  AOA. 

Clear  [h]; 

Clear  [theta]; 

Clear  [equations]; 

Clear  [coeffs]; 

Clear  [sol];  %  Clear  all  variables 

Clear  [Theta]; 

Clear  [H]; 

Clear  [aeff]; 


Youngs  =  Input  [’’Young’s  Modulus  =  ”]; 
G  =  Input  [’’Modulus  of  Rigidity  =  ”]; 
pois  =  Input  [’’Poisson’s  Ratio  =  ”]; 

El  =  Youngs/G; 

comp  =  Input  [”1)  Compressible 

Analysis,  2)  Incompressible  Analysis”]; 

alphain  =  Input  [’’Angle  of 
Attack(deg)  =  ”]; 
alphaO  =  alphain  *  3.14159/180; 
alphabar  =  alphaO; 
alphadbar  =  alphaO; 

c  =  Input  [’’Chord  (ft)  =  ”]; 
b  =  Input  [’’Wing  Span  (ft)  =  ”]; 

1  =  b/2; 

e  =  Input  [’’Elastic  Offset  (ft)  =  ”]; 
ebar  =  e/c; 

t  =  Input  [’’Wing  Thickness(ft)  =  ”]; 
tbar  =  t/c; 

Sweep  =  Input  [’’Wing 
Sweepback(deg)  =  ”]; 
sweep  =  Sweep  *  3.14159/180; 

AR  =  b/c; 

ARl  =  AR/2; 

V  =  Input  [”Airspeed(ft/s)  =  ”]; 
roe  =  Input[”Density(slugs/ft^)  =  ”]; 

T  =  Input  [’’Temperature  (deg  F)  =  ”]; 
sound  =  (1.4*1716*(T  +  460))'5; 

M  =  V/sound; 


%  Input  Young’s  modulus 
%  Input  shear  modulus 
%  Input  Poison’s  ratio 
%  Calculate  E/G’ 
%  Choose  flow  regime 

%  Input  AOA  in  deg 

%  Degree  to  radians 
%  Angle  of  attack  (1st  normalization) 
%  Angle  of  attack  (2nd  normalization) 

%  Input  chord 
%  Input  wingspan 
%  Calculate  semispan 
%  Input  elastic  offset 
%  Calculate  elasic  offset 
%  Input  wing  thickness 
%  Calculate  wing  thickness 

%  Input  sweep 

%  Degree  to  radians 
%  Calculate  aspect  ratio 
%  Calculate  semi-aspect  ratio 
%  Input  velocity 
%  Input  density 
%  Input  temp 
%  Calculate  speed  of  sound 
%  Calculate  Mach  number 
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qn  =  .5  *  roe  *  *  Cos  [sweep] 

Qn  =  qn  *  12  *  (1  -  pois^)  * 

ARl/ (Youngs  *  tbar^); 

aincomp  =  6. 28/(1  +  (4/AR)  Cos  [sweep]); 
acomp  =  (aincomp*Cos[sweep])/((l  - 
M^*Cos[sweep]^  +  (aincomp*Cos[sweep]/ 
(3.14159*AR))2)-5  + 
aincomp*  Cos  [sweep]  /  (3. 14159*  AR) ) ; 

cmac  =  0; 

If[comp  ==  2,  aO  =  aincomp,  aO  =  acomp]; 

ml  =  (tbar)^/(12  (1  -  (pois)^)  (ARl)^); 

m2  =  (tbar)^/(2  (1  +  pois)); 

m3  =  ml  aO  ARl; 

m4  =  ml  aO  Tan  [sweep]; 

m5  =  ml  aO  alphabar  ARl  Qn; 

m6  =  ((tbar)^/(ARl  (1  -  (pois)^)))*ebar  aO; 

m7  =  m6  Tan  [sweep] /ARl; 

m8  =  m6  alphadbar  Qn; 

m9  =  ((tbar)^/((l  -  (pois)^)  ARl))*cmac  Qn; 

ml2  =  m2/ml; 
ml3  =  m3/ml; 
ml4  =  m4/ml; 
ml5  =  m5/ml; 
ml6  =  m6/ml; 
ml7  =  m7/ml; 
ml8  =  m8/ml; 
ml9  =  m9/ml; 

al  =  -m4  El  Qn; 
a2  =  m3  El  Qn; 
a3  =  -ml3  Qn; 
a4  =  ml4  Qn; 
bl  =  ml5; 


%  Calculate  dynamic  pressure 
%  Normalize  dynamic  pressure 

%  Incompressible  lift  curve  slope 
%  Compressible  lift  curve  slope 

%  Moment  coefficient 
%  Choose  ffow  regime 

%  Define  coefficients  from  governing  equations 


%  Define  coefficients  from  governing  equations 


%  Define  coefficients  from  governing  equations 
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a5  =  -(4  ml2  -  m6  El  Qn)/(m2  El  +  1); 
a6  =  -ml6  Qn; 
a7  =  ml7  Qn; 

a8  =  -(m7/(m2  El  +  1))  El  Qn; 
b2  =  mis  +  nil9; 

cl  =  ml  El; 

c2  =  -ml  m4  El^  Qn; 

c3  =  1; 

c4  =  ml  m3  El^  Qn; 

c5  =  ml  El  (m2  El  -|-  1); 

c6  =  ((m2  El  -|-  1)^  -|-  ml  m6  El^  Qn); 

c7  =  -ml  m7  El^  Qn; 

dl  =  m3  El  Qn; 

d2  =  -m4  El  Qn; 

b3  =  -m5  El; 

d3  =  (m2  El  +  1);; 

d4  =  m6  El  Qn; 

d5  =  -m7  El  Qn; 

b4  =  -(m8  -|-  m9)  El; 

d6  =  -m4  El  Qn; 

d7  =  m3  El  Qn; 

d8  =  (m2  El  +  1); 

d9  =  -(4  ml2  -  m6  El  Qn); 

dlO  =  -m7  El  Qn; 

sol  =  DSolve[ 

{  h””[y]  +  al  h”’[y]  +  a2  theta”  [y] 

-|-  a3  theta[y]  -|-  a4  h’[y]  ==  bl, 
theta”  ”  [y]  -|-  a5  theta”  [y] 

-|-  a6  theta[y]  -|-  a7  h’[y]  -f  a8  h”’[y]  == 
b2},  {h[y],  theta[y]},  y]; 


%  Dehne  coefficients  from  governing  equations 


%  Dehne  coefficients  from  governing  equations 


%  Solve  governing  equations 


%  Solve  governing  equations 
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h[y_]  =  First [h[y]  /.  sol]; 
theta[y_]  =  First  [theta  [y]  /.  sol]; 
equations  =  {h[0]  ==  0, 
theta[0]  ==  0, 
cl  h’”[0]  +  c2  h”[0] 

+  c3  h’[0]  +  c4  theta’ [0]  ==  0, 
c5  theta” ’[0]  +  c6  theta’ [0] 

+  c7  h”[0]  ==  0, 

h”[l]  +  dl  theta[l]  +  d2  h’[l]  ==  b3, 
d3  theta”  [1]  +  d4  theta[l]  , 

+  d5  h’[l]  ==  b4 

h’”[l]  +  d6  h”[l]  +  d7  theta’ [1]  ==  0, 
d8  theta” ’[1]  +  d9  theta’ [1] 

+  dlO  h”[l]  ==  0}  ; 

%  Input  boundary  conditions 

coeffs  =  Solve  [equations,  {C[l],  C[2],  C[3], 
C[4],  C[5],  C[6],  C[7],  C[8]}]; 

%  Solve  for  unknowns 

Theta[y_]  =  Chop  [ComplexExpand  [Simplify 
[First[theta[y]  /.  coeffs]]]]; 

%  Solve  for  twist 

H[y_]  =  Chop  [ComplexExpand  [Simplify 
[First [h[y]  /.  coeffs]]]]; 
aeff[y_]  =  1  +  (Theta[y]  - 

%  Solve  for  bending 

(Tan[sweep]/AR1)  D[H[y],  {y,  l}])/alpha0; 

%  Solve  for  effective  AOA 

plotA  =  Table[{l*y,  Re[H[y]*c]}, 

{y,  0,  1,  .01}]; 

%  Eliminate  imaginary  term 

plotB  =  Table[{l*y,  Re[Theta[y]]}, 

{y,  0,  1,  .01}]; 

%  Eliminate  imaginary  term 

plotC  =  Table[{l*y,  Re[aeff[y]]*alphain}, 

{y,  0,  1,  .01}]; 

%  Eliminate  imaginary  term 

ListPlot  [plotA,  Plot  Joined  -i  True, 

AxesLabel {’’Non-Dimensional  Semispan”, 
’’Vertical  Displacement  (ft)”}]; 

%  Plot  bending 

ListPlot  [plotB,  Plot  Joined  -i  True, 

AxesLabel {’’Non-Dimensional  Semispan”, 
’’Twist  (deg)”}]; 

%  Plot  twist 

ListPlot  [plotC,  Plot  Joined  -i  True, 

AxesLabel {’’Non-Dimensional  Semispan”, 

%  Plot  effective  AOA 

’’Effective  Angle  of  Attack  (deg)”}]; 
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Table  A. 6:  Unsteady  aerodynamic  coefficients. 

Lhh  =  2  (^1  +  2^^  — 

Lhe  =  2a  +  2A  -  4^  (^i  -  a)  +  I  [(l  +  2f  +  4  (I  -  a))]  i 
Lhh'  =  -4p^  tan  A  -  ^  tan  A  (l  +  2f )  * 

Lhe'  =  (I  -  a)  tan  A  -  ^  tan  A  [a  -  |  -  a)  G]  i 

Mgh  =  4a  +  8^  (|  +  ~  Sf"  (|  +  ®)  * 

Mee  =  l+  4a2  +  4^  Q  +  a)  -  8f  ( A  -  a^)  -  f  [4  -  a  -  2f  (l  +  a)  -  AF  (A  _  a^)]  i 
^eh'  =  8p^  tanA  {\  +  a)-^  tan  A  [a  +  2f  (4  +  a)]  i 
^oe'  =  p^tanA[l-16F(A-a2)]  _  ^  tan  A  [A  +  -  2f  (A  _  a^)]  i 


